Bold Diagrammatic Monte Carlo: When Sign Problem is Welcome 
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We introduce a Monte Carlo scheme for sampling bold-line diagrammatic series specifying an 
unknown function in terms of itself. The range of convergence of this bold (-line) diagrammatic 
Monte Carlo (BMC) is significantly broader than that of a simple iterative scheme for solving 
integral equations. With BMC technique, a moderate "sign problem" turns out to be an advantage 
in terms of the convergence of the process. For an illustrative example, we solve one-particle s- 
scattering problem. As an important application, we obtain T-matrix for a fermipolaron (one 
spin-down particle interacting with the spin- up fermionic sea). 
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Diagrammatic Monte Carlo (DiagMC) [ll] is a tech- 
nique that allows one to simulate quantities specified in 
terms of convergent diagrammatic sums, i.e., sums of in- 
tegrals with integrands represented by a diagrammatic 
structure. Formally, it is a set of generic prescriptions for 
organizing a systematic-crror-free Metropolis-type pro- 
cess that samples the series/sum without explicitly tak- 
ing the integrals over the internal variables in each par- 
ticular term. In addition to the natural requirement of 
convergence, the diagrammatic sums should be either es- 
sentially finite (have only a few leading terms) or positive 
definite; otherwise the sign problem suppresses the effi- 
ciency of the numeric procedure. One of the key tools 
in the analytical diagrammatic techniques is the trick of 
bold lines [3| that allows one to (partially or completely) 
sum the scries even if it is formally divergent. The bold- 
line trick looks also very attractive for the sign-indefinite 
series since it can substantially reduce the number of 
leading diagrams, and thus alleviate the sign problem. 

In this Letter, we explore the possibility of employing 
the bold-line trick in the DiagMC approach. We propose 
a scheme which we call bold(-linc) diagrammatic Monte 
Carlo (BMC). In essence, BMC is a generalized iterative 
scheme in which the iteration protocol depends on the 
number of iteration steps, or, equivalently, in which the 
next iteration is a function of not only its immediate an- 
cestor, but of the (properly weighted) whole list of earlier 
iterations. The crucial difference between BMC and a 
nai've iteration protocol — when one simply uses DiagMC 
to perform an integration for a given iteration step — is 
that the convergence of BMC has essentially broader pa- 
rameter range. We present general arguments why this is 
the case and perform an illustrative simulation for one- 
particle s-scattering problem. Despite its formal simplic- 
ity, the problem contains all the ingredients one can meet 
in a general case: the pcrturbativc series diverges if the 
scattering potential is strong enough, and — in the case of 
a repulsive potential — the series is not positive definite. 
The simplifications which we have here are mainly quan- 
titative rather than qualitative. The bold-line trick re- 
duces the infinite perturbative series to just two terms. In 
the case of a strong attractive potential, two more terms 



appear in the right-hand side to secure the convergence. 
Incidentally, with these extra terms a sign problem arises 
for strong attractive potential as well. 

Before turning to the realistic simulation, we start with 
discussing a generalized iterative procedure which is most 
close to the actual BMC scheme. To analyze the conver- 
gence issues, we can confine ourselves with a linearized 
problem and write 



I/) = |6) + A|/), 



(1) 



where | / ) is some unknown vector, | 6 ) is a known vector, 
and A is a linear operator. Expanding all the vectors in 
terms of the operator A eigenvector basis {| )}, we get 



.h = h + fi 



(2) 



where A| 



= ai\h), \b) = Eih\h); l/) = 
) . The vector equation thus decouples into a 
set of independent equations for each f . From now on the 
label ^ can be omitted. From the convergence point of 
view, it is advantageous to work with iterative schemes 
that involve — at least at the theoretical level — an aver- 
aging of the iterations. Let /„ be the n-th generator for 
the (n -|- l)-st iteration 
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The quantity /„ is supposed to be a function of all /j 's 
with J < n. As a characteristic example we choose 



fn 



(4) 



where a > —1 is a fixed parameter of the scheme. We can 
exclude /'s and explicitly relate f(n+i) to /„ to see that 
for the deviation i5n = / — /n of the n-th generator from 
the exact solution / — 6/(1 — a) the following recursive 
relation takes place 



^(n+l) 



= Sn 



(l + a)(a-l) 



It implies the asymptotic behavior 
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Hence, the condition of convergence is 
Rea^ < 1 . 



(7) 



Here we restore the subscript ^ to emphasize that condi- 
tion ([7]) has to be satisfied for all the eigenvalues of the 
matrix A. We see that the value of a does not determine 
the fact of convergence, but does effect the asymptotic 
rate of convergence — the larger is a, the higher the rate. 
It is important that the convergence does not depend on 
the imaginary parts of a^'s. Finally, negative real parts 
of are desirable for convergence: the larger is the ab- 
solute value of the negative real part, the better. [Note 
that the plain iterative method (/„ = /„) converges only 
when \a^\ < 1.] 

If condition ([7]) is not met. one can use an equation 
equivalent to ([l]), but with convergent iterative proce- 
dure. For the s-scattering problem, the matrix A is Her- 
mitian and thus all its eigenvalues are real. In this case, 
rewriting Eq. ([T]) as 

/) = |fe)+A|/) + AA(|/)-|6)-A|/)) (8) 

and fine-tuning the value of the constant A, one can 
render the iterative process converging. Indeed, the 
new equation has the same form as the original one, 
up to replacements \b) — > |^) — AA|&) and A — > 
(1 -f A)A — AA^. Correspondingly, condition ([7]) for the 
new matrix will be met if the original eigenvalues satisfy 
the inequality 



(1 4- \)a^ - Aa| < 1 



(9) 



As is easily checked, condition ^ is met provided A G 
(Ai, A2), where 



\ = min{aj 

aj>l 



An ^ max jflf I . (10) 



Hence, if the eigenvalues are real and separated from 
unity by a finite gap, there exists a value of A at which 
convergence is guaranteed. Incidentally, the problem ([T|) 
can always be re-formulated in such a way that the new 
matrix is Hermitian: 

/) = (1-At)|6) + (A + At-AtA)|/). (11) 

The s-scattering problem can be formulated (see, e.g., 
Ref. 3; for simplicity, we work with a spherically sym- 
metric potential) as Eq. ^ with |/) = f{q) being 
the zero-energy scattering wave function in the momen- 
tum representation. In this case |&) = ~u{q), where 
u{q) = U{q)/2-K, and U{q) is the Fourier transform of 
the scattering potential; the Plank's constant and parti- 
cle mass are set equal to unity. The operator A here is 
the integral operator 



A/ = - - 



- I dxf ti(lq-qil) /(gi)rfgi , (12 

J -I Jo 



where |q — qi| = \/ + — 2qqix; Eq. ^ then reads 
f = -u + XAu + (1 + A)A/ - AAV . (13) 



The potential we use in simulations is a flat spherical 
well/bump defined as J7(r) = C/q at r < 1 and zero oth- 
erwise. For this potential. 



HQ) = — 5- (smq - gcosq) 
q-^ 



(14) 



Monte Carlo procedure. Here wc discuss how generic Di- 
agMC rules are used to calculate f{q) self-consistently. 
For brevity, let us index the four terms in the right-rand- 
side of Eq. as terms A, B, C, and V. Correspond- 
ingly, the "configuration space" of the problem is defined 
by the term index and all continuous variables associated 
with it. The goal of the Monte Carlo procedure is to per- 
form stochastic summation over this configuration space. 
The contribution of each state to the answer is charac- 
terized by the weight with the sign which in our case is 
given by the product of u- and /-functions; for example, 
the weight and sign of the (B, q, q' , x) state is determined 
by the modulus and sign of u{\q — q'\)f{q')- 

The standard Metropolis-type protocol consists of up- 
dates which change the current configuration state fol- 
lowed by measurements which evaluate contributions of 
the current state to the answer. The updating scheme de- 
scribed below generates states with probabilities propor- 
tional to their weight. In this case, the Monte Carlo esti- 
mator for f{q) is given by the state sign. The statistics of 
±1 contributions is accumulated into the /((3')-histogram 
with bins covering the positive-g axis. Apart from repre- 
senting the final result of the simulation, the histogram is 
used self-consistently to generate random variables from 
the probability density \f{q)\ and to determine the sign 
of B and V states. 

A straightforward accumulation of data into the his- 
togram corresponds to a = 1 in Eq. ([4]). However, large 
values of a result in faster convergence, see Eq. 1^. Nu- 
merically, the limit of a ^ 00 is implemented by simply 
erasing "old" histogram data. 

The simplest updating scheme contains three pairs of 
complementary updates [A ^ B,B ^ A], [A ^ C,C ^ 
A], [C — > — > C] which change the term index, and 
one self-complementary update A ^ A changing the 
variable q. 

A ^ A update. A new value for the variable q in state 
A is generating from the normalized probability density 



p{q) = \u{q)\/Iu 



\u{q)\dq , 



(15) 



The acceptance ratio for the A A update is unity. 

A ^ B update. First, we select the value of x from 
the uniform probability density on the [—1,1] interval. 
Next, we select the value q' from the histogram based 
probability distribution \f{q')\. The acceptance ratio for 
this update is 



Ra^b 



2|1 + A|// 

T^PAB 



u{q - q') 



u{q) 



(16) 



where pab is the probability to apply the A ^ B update 
while in state A (we do not mention the probability of 
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applying an update to the current configuration if it is 
unity; in this scheme pbj\ = !)• 



If = 



\fi<l)\dq 



(17) 



is the normalization integral proportional to the sum of 
absolute values of all histogram contributions (its proper 
normalization is discussed below). 

B A update. Here we propose to change the term 
index back to A] the acceptance ratio for this move is 
simply the inverse of Rji^js. 

A C and C A updates. Formally, this pair of 
updates is identical in implementation to the previous 
one with the only difference that the value of q' in the 
A —> C move is generated from the probability density 
|u((7')|. Correspondingly, the acceptance ratio is based 
on the lu integral: 
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2|A| luPcA 



u{q - q') 



u{q) 



(18) 



Rc^A T^PAC 

T) and V ^ C updates are an exact copy of the 
A ^ B and B ^ A pair in terms of how new variables 
are proposed and removed. The acceptance ratio is 
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T^PCV 



u{q' 



u{q') 



(19) 



The above set of updates is ergodic, i.e. it sam- 
ples the entire configuration space. In the practical 
implementation of the algorithm we used PAA ~ 0.2, 
PAB = PAC = 0.4, and pcA = Pct> = 0.5. To complete 
the description we have to explain how to find the nor- 
malization integral // in Eq. (HI]). Let Zmc be the total 
number of Monte Carlo states in the simulation and Zy^ 
be the number of ^-states. In the statistical limit. 



Za 

Zmc 



In 

I 



(20) 



I^If + 



|1 + A| 



\u{q-q')f{q')\dxdqdq' 
+ V/ - 9>il')\dxdqdq' (21) 
+ 1^/ Hq-q')u{q' -q")f{q")\dxdx'dqdq'dq" . 

is the auxiliary "global partition function" which drops 
out from all final answers. If Hg is the sum of all con- 
tributions to the s-th bin of the histogram then in the 
statistical limit, 



Jk 

Zmc 



= r 



f{q)dq. 



(22) 



If we now write the normalization integral as a sum over 
the histogram (in the limit of infinitesimally small bin 
size the relation is exact) 



fo{q) dq 



{Z/Zmc) E I^^I ' (23) 



and use Eq. ([20]) to eliminate I /Zmc we finally arrive at 



Za 



(24) 



Similarly, the physical result for the scattering wave func- 
tion is given by 



f{qs) §^i^ ■ 



(25) 



The s-wave scattering length can be obtained in two 
ways: from the q —> limit, a ~ —f{q ~ 0), and as a 
histogram sum (the last procedure gives better accuracy 
since it is based on all Monte Carlo data and thus is not 
susceptible to the noise in a particular bin) 



a = u{0) + - 

TT 



u{q)f{q)dq ^ w(0)- 



2Iu 



^ u{qs)Hs 



(26) 

We have tested our BMC scheme against the analyti- 
cal answer for the scattering length in different regimes 
which included strong repulsive and attractive potentials 
outside of the convergence limits for the standard iter- 
ative scheme. For example, one can easily get results 
for a with four-digit (or higher) accuracy for the repul- 
sive potential Uq = 10; a straightforward summation of 
the series expansion for large positive values of Uq would 
be impossible because of the divergence and/or the sign 
problem. Series divergence will also prevent one from 
going across the resonance and getting results for poten- 
tials with bound states. In Fig. [T] we present data for 
the scattering wave function obtained for Uo = —3, i.e. 
for the potential well with the bound state. In this sim- 
ulation a was obtained with the 4-digit accuracy. For 
negative values of Uq < —10 we found that good initial 
conditions, e.g. results of the previous run for smaller 
\Uo\, are important for convergence which was very slow 
and required that Awl. In view of Eq. pU]) . we con- 
clude that in this parameter range we face the problem 
of having matrix A eigenvalues being too close to unity. 

Fermipolaron T-matrix. The fermipolaron is a spin- 
down particle interacting with the sea of spin-up 
fermions. Of special interest is the case when the spin-up 
sea is an ideal Fermi gas while the interaction between 
spin-up and spin-down particles is short-ranged but res- 
onantly strong. In this regime — relevant to the notorious 
problem of BCS-BEC crossover in the limit of extreme 
population imbalance between the two fermionic compo- 
nents ^ — there is, in particular, a critical point in the 
interaction strength when the groundstate of the polaron 
becomes a bound spin-zero state (molecule). The fer- 
mipolaron problem allows an unbiased numeric solution 
by DiagMC. The relevant diagrams are constructed out 
of the three types of propagators (in the imaginary-time- 
momentum representation): (i) spin-up Green's func- 
tions, (ii) spin-down vacuum propagators, and (iii) the T- 
matrix of the pair interaction between spin-up and spin- 
down particles Q. While simple analytic expressions for 
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FIG. 1: (Color online). Scattering wave function at zero en- 
ergy (solid line) and scattering potential (dashed line) for the 
attractive potential well with one bound state {Uo = —3). 



FIG. 2: The diagrammatic equation for the T-matrix (heavy 
dashed line) in terms of the vacuum T-matrix (light dashed 
line), spin-down vacuum propagator (straight solid line), and 
truncated (to the momenta less than Fermi momentum) spin- 
up propagator (solid arc). 

type (i) and (ii) propagators are available, the T-matrix 
has to be tabulated numerically; this tabulation repre- 
sents the performance bottleneck for the whole scheme. 
Noting that the equation defining the T-matrix through 
itself (sec Fig. [2]) is analogous to that for the scattering 
amplitude, one can directly apply the above-described 
BMC procedure for obtaining the T-matrix. We have 
successfully done that, which ultimately allowed us to 
solve the fermipolaron problem Q . 

Conclusions and outlook. We have found a numeric 
counterpart to the bold-line trick of diagrammatic tech- 
nique. The resulting scheme simulates unknown func- 
tions specified by diagrammatic series in terms of them- 
selves. We illustrated our approach by solving the s- 
scattering problem in strong repulsive and attractive po- 



tentials. We introduced tools to secure convergence of 
the process. With these tools we were able to solve the 
s-scattering problem even in an attractive potential with 
a bound state — an essentially non-groundstate problem. 

The standard many-body diagrammatic technique 
deals with three functions that are expressed through 
each other: Green's function, self-energy, and the four- 
point vertex. The generalization of the scheme to this 
case is theoretically straightforward, since formally one 
can think of all these functions as different components 
of the vector |/). The three practical questions that are 
immediately seen — in the order of their importance — are: 
(i) regularization of bold-line series, (ii) convergence of 
the scheme, (iii) optimal data structure. Formally, the 
convergence of the bold-line series may depend on the 
summation order and in certain cases be achievable at 
the expense of controllable systematic error. One may 
keep the expansion-order, m, of irreducible diagrams fi- 
nite and extrapolate results to the to — > oo limit, work 
with a finite-size system and extrapolate to the thermo- 
dynamic limit, introduce constraints on continuous vari- 
ables not allowing them to be either very small or very 
large and apply renormalization techniques (ultra-violet 
divergences would be a typical example), etc. The con- 
vergence of the scheme may be achieved by the tools 
described in this paper. If the initial approximation to 
unknown functions is close enough to the exact answer — 
which will be the case if one starts with an almost ideal 
system and moves to a strongly interacting regimes by 
small steps in the interaction constant, then one may 
rely on linearization for constructing the correcting part 
of the right-hand side operator, using prescriptions of 
Eq. ([8]). If unknown functions depend on many con- 
tinuous variables, histograms may well be not the best 
method. Instead one may use a variable-step meshes 
and, correspondingly, re-weighing techniques for collect- 
ing statistics. Another option is to approximate unknown 
functions with analytic expressions and permanently op- 
timize their parameters to best fit coarse-grained his- 
togram sampling coming from the BMC process. 
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